*变量描述
local varlist "eco7 indus eco13 eco1 city15 env65  env59 env93 yr_sch_sec  "
tabstat env38  wbgi_rqe  `varlist', s(n mean median sd min max) f(%9.4f)  c(s)

* 散点图
	* 区分1998年签订国和未签订国
twoway   ( scatter env38 wbgi_rqe if sign !=1 | countryname!=, msymbol(T)  msize(0.5)) ///
		  ( qfit env38 wbgi_rqe if sign !=1 ) ///
		( scatter env38 wbgi_rqe if sign ==1 ,msymbol(o) msize(0.5)) ///
	    (qfit env38 wbgi_rqe if sign ==1  ) 
	* 区分2002年时间节点
twoway   ( scatter env38 wbgi_rqe if year<2003 , msymbol(T) msize(0.5)) ///
		  ( lfit env38 wbgi_rqe if year<2003 , color(gray)) ///
		( scatter env38 wbgi_rqe if year>2002 ,msymbol(o) msize(0.5)) ///
	    (lfit env38 wbgi_rqe if year>2002 )

twoway   ( scatter env38 wbgi_rqe if year<2003 , msymbol(T) msize(0.5)) ///
		  ( qfit env38 wbgi_rqe if year<2003 , color(gray)) ///
		( scatter env38 wbgi_rqe if year>2002 ,msymbol(o) msize(0.5)) ///
	    (qfit env38 wbgi_rqe if year>2002 )
	
 

* 按2018年世界银行收入分类
			
			gen incomelev=.
			replace incomelev=1 if (income=="H" & year==2018)
			gen incomelev1 = incomelev if year ==2018
			bysort countrycode: egen incomelev2=min(incomelev1)
			
    keep if (year >1994 & year<2015)
	set scheme  s2mono
twoway   ( scatter env38 eco7 if incomelev2!=1 , msymbol(o) msize(0.5)) ///
		  ( lfit env38 eco7 if incomelev2!=1, color(red))
		  
		
twoway   ( scatter env38 eco7 if incomelev2==1 , msymbol(o) msize(0.5)) ///
		  ( lfit env38 eco7 if incomelev2==1, color(red) ylabel(0(0.5)2) )

twoway   ( scatter env38 wbgi_rqe if incomelev2!=1 , msymbol(o) msize(0.5)) ///
		  ( lfit env38 wbgi_rqe if incomelev2!=1, color(red))
		  
		
twoway   ( scatter env38 wbgi_rqe if incomelev2==1 , msymbol(o) msize(0.5)) ///
		  ( lfit env38 wbgi_rqe if incomelev2==1, color(red))
		  
reg env38 wbgi_rqe eco1 eco13 indus city15 env65  env59 env93 yr_sch_sec
vif		  
		xtset id year  
	local varlist "  eco1  city15 env65  env59 env93 yr_sch_sec "  
	reg env38 wbgi_rqe if incomelev2!=1
   est store OLS2
   reg env38 wbgi_rqe if incomelev2==1
   est store OLS3
   xtreg env38 wbgi_rqe  `varlist' if incomelev2!=1,fe // 113 96-14
   est store fe1
   xtreg env38 wbgi_rqe   `varlist' if incomelev2==1,fe  // 113 03-14
   est store re
   esttab  OLS2 OLS3 fe1 re ,  label ///
    compress nogap  s(N r2  F)  ///
          star(* 0.1 ** 0.05 *** 0.01)  b(%6.4f) 	
		  

*3.OLS,面板回归基本结果
* 静态面板回归
gen eco7w =eco7/10000
preserve
keep if  (year >1994 & year <2015)
xtbalance,range(1995,2014)
replace eco7w = L.eco7w if eco7w==.
xtunitroot llc eco7w
restore
local varlist "eco7w eco13 eco1 city15 env65 env24 env59 env93 yr_sch_sec hf_efiscore"
drop id 
egen id = group(countrycode)
xtset id year
gen rqe2 = wbgi_rqe^2



* 系数扩大100倍
gen envx = env38*100
drop env38
rename envx env38

gen envx = env38/100
drop env38
rename envx env38


* Kyoto区分回归
local varlist "eco7w indus eco1 eco13  city15 env65   env59 env93 yr_sch_sec "
reg env38 wbgi_rqe if sign==1
est store OLS1
reg env38 wbgi_rqe if sign !=1
est store OLS2
xtreg env38 wbgi_rqe `varlist'  if sign==1,fe // 113 03-14
est store fe1
xtreg env38 wbgi_rqe `varlist' i.year if sign==1,fe // 113 03-14
est store fe2
xtreg env38 wbgi_rqe   `varlist' if sign !=1,fe // 113 03-14
est store fe3
xtreg env38 wbgi_rqe   `varlist' i.year if sign !=1,fe // 113 03-14
est store fe4
esttab OLS1 OLS2  fe1  fe2 fe3 fe4, ///
 compress nogap  s(N r2  F)  ///
          star(* 0.1 ** 0.05 *** 0.01)   b(%6.4f) 	
	
*ratification区分回归	
drop ratif
gen ratif = .
replace ratif = 1 if year > ratification
local varlist "eco7w indus eco1 eco13  city15 env65  env59 env93 yr_sch_sec "
reg env38 wbgi_rqe if (ratif ==1  & countrycode!="USA" & countrycode != "CAN")
est store a1
reg env38 wbgi_rqe if (ratif !=1 & countrycode!="USA" & countrycode != "CAN")
est store a2
reg env38 wbgi_rqe `varlist' if (ratif ==1 & countrycode!="USA" & countrycode != "CAN")
est store x1 
reg env38 wbgi_rqe `varlist' if (ratif !=1 & countrycode!="USA" & countrycode != "CAN")
est store x2
esttab a1 a2 x1 x2, ///
 compress nogap  s(N r2  F)  ///
          star(* 0.1 ** 0.05 *** 0.01)
		  

		  
* DID
drop if year <1995
bysort year: egen env38m=mean(env38)	if sign == 1
bysort year: egen env38n=mean(env38) if sign!=1

label var env38m "treatment group"
label var env38n "control group"



drop did
gen did = .
replace did = 1 if year >2002
replace did = 0 if year <2003
replace sign = 0 if sign !=1
local varlist "eco7w indus eco13 eco1 city15 env65  env59 env93 yr_sch_sec "
xtreg env38 wbgi_rqe did  sign#did ,fe   // 240 95-14
est store did1
xtreg env38 wbgi_rqe did  sign#did  i.year,fe    // 240 95-14
est store did2
xtreg env38 wbgi_rqe did  sign#did `varlist',fe   // 240 95-14
est store did3
xtreg env38 wbgi_rqe did  sign#did `varlist' i.year,fe    // 240 95-14
est store did4

esttab did1  did2 did3 did4 ,compress nogap  s(N r2  F ) ///
          star(* 0.1 ** 0.05 *** 0.01)   b(%6.4f)
		

twoway   ( scatter env38 wbgi_gee if sign !=1 , color(gray) msize(0.5)) ///
		  ( lfit env38 wbgi_gee if sign !=1 , color(bulue)) ///
		( scatter env38 wbgi_gee if sign ==1 , color(blue) msize(0.5)) ///
	    (lfit env38 wbgi_gee if sign ==1 , color(blue)) 
		

twoway   ( scatter env38 lpat if sign !=1 , color(gray) msize(0.5)) ///
		  ( lfit env38 lpat if sign !=1 , color(bulue)) ///
		( scatter env38 lpat if sign ==1 , color(blue) msize(0.5)) ///
	    (lfit env38 lpat if sign ==1 , color(blue)) 
		
		
		
		
* regtype分类回归

			*按2008年的regime type分类
			gen chga_hinst1 = chga_hinst if year ==2008
			bysort countrycode: egen chga_hinst2=min(chga_hinst1)	
			drop if year<1995
			
			xtset id year
local varlist "eco7w eco1  indus eco13  city15 env65  env59 env93 yr_sch_sec "
xtreg env38 wbgi_rqe `varlist' if chga_hinst2 == 0,fe
est store parliamentary
xtreg env38 wbgi_rqe `varlist' if  chga_hinst2  == 1,fe
est store mixed
xtreg env38 wbgi_rqe `varlist' if chga_hinst2  == 2,fe
est store presidential
xtreg env38 wbgi_rqe `varlist' if (chga_hinst2  == 3 | chga_hinst2==4 | ///
									chga_hinst2==5),fe
est store dictatorship

esttab parliamentary mixed presidential dictatorship,compress nogap  s(N r2  F ) ///
          star(* 0.1 ** 0.05 *** 0.01)   b(%6.4f)
		  
		  
			*筛选所有区间regime type没用变化的
			gen chga_hinstp = .
			bysort countrycode: egen hinstmax=max(chga_hinst)
			bysort countrycode: egen hinstmin=min(chga_hinst)
			bysort countrycode: replace chga_hinstp = hinstmax if hinstmax==hinstmin
			drop if year<1995
			
			
			*筛选从1998-2008年regime type没用变化的国家
			sort id year
			drop hinstmax
			gen hinstmax=100
			bysort countrycode: replace hinstmax=max(chga_hinst[4],chga_hinst[5],chga_hinst[6],chga_hinst[7],chga_hinst[8],chga_hinst[9],chga_hinst[10],chga_hinst[11],chga_hinst[12],chga_hinst[13],chga_hinst[14])
			sort id year
			drop hinstmin
			gen hinstmin=-100
			bysort countrycode: replace hinstmin=min(chga_hinst[4],chga_hinst[5],chga_hinst[6],chga_hinst[7],chga_hinst[8],chga_hinst[9],chga_hinst[10],chga_hinst[11],chga_hinst[12],chga_hinst[13],chga_hinst[14])
			
			drop chga_hinstp
			gen chga_hinstp = .
			bysort countrycode: replace chga_hinstp = hinstmax if hinstmax==hinstmin
			drop if year<1995
			
			*regime type作图
			set scheme s2mono
	twoway   ( scatter env38 wbgi_rqe if ///
		(chga_hinstp == 0 | chga_hinstp==1 |chga_hinstp==2) , ///
		msymbol(O) msize(0.5)) 	  ( qfit env38 wbgi_rqe if ///
		(chga_hinstp == 0 | chga_hinstp==1 |chga_hinstp==2) ) ///
 ( scatter env38 wbgi_rqe if (chga_hinstp  == 3 | chga_hinstp==4 | ///
									chga_hinstp==5) , msymbol(T) msize(0.5)) ///
	    (qfit env38 wbgi_rqe if (chga_hinstp  == 3 | chga_hinstp==4 | ///
									chga_hinstp==5) , ) 

twoway ( scatter env38 wbgi_gee if dem ==1 ,msize(0.2) msymbol(O) color(maroon) )(lfit env38 wbgi_gee if dem==1, color(red)) ( scatter env38 wbgi_gee if dem !=1 ,msize(0.2) msymbol(T) color(gray))(lfit env38 wbgi_gee if dem !=1, color(blue))			
												
			*regime type回归						
		  local varlist "eco7w eco1  indus eco13  city15 env65  env59 env93 yr_sch_sec "
xtreg env38 wbgi_rqe `varlist' if ///
		(chga_hinstp == 0 | chga_hinstp==1 |chga_hinstp==2),fe
est store m1
xtreg env38 wbgi_rqe `varlist' i.year if ///
		(chga_hinstp == 0 | chga_hinstp==1 |chga_hinstp==2),fe
est store m2
xtreg env38 wbgi_rqe `varlist' if (chga_hinstp  == 3 | chga_hinst2==4 | ///
									chga_hinst2==5),fe
est store m3
xtreg env38 wbgi_rqe `varlist' i.year if (chga_hinstp  == 3 | chga_hinst2==4 | ///
									chga_hinst2==5),fe
est store m4

esttab m1 m2 m3 m4,compress nogap  s(N r2  F ) ///
          star(* 0.1 ** 0.05 *** 0.01)   b(%6.4f)
	
	*加入二次项的回归
	gen gee2=wbgi_rqe^2
  local varlist "eco7w eco1  indus eco13  city15 env65  env59 env93 yr_sch_sec "
xtreg env38 wbgi_rqe gee2 if (chga_hinstp == 0 | chga_hinstp==1 |chga_hinstp==2),fe
est store ad1
xtreg env38 wbgi_rqe gee2 if (chga_hinstp== 3 | chga_hinstp==4 |chga_hinstp==5),fe 
est store ad2
xtreg env38 wbgi_rqe gee2 `varlist' if ///
		(chga_hinstp == 0 | chga_hinstp==1 |chga_hinstp==2),fe
est store m1
xtreg env38 wbgi_rqe gee2 `varlist' i.year if ///
		(chga_hinstp == 0 | chga_hinstp==1 |chga_hinstp==2),fe
est store m2
xtreg env38 wbgi_rqe gee2  `varlist' if (chga_hinstp  == 3 | chga_hinstp==4 | ///
									chga_hinstp==5),fe
est store m3
xtreg env38 wbgi_rqe gee2 `varlist' i.year if (chga_hinstp  == 3 | chga_hinstp==4 | ///
									chga_hinstp==5),fe
est store m4

esttab ad1 ad2 m1 m2 m3 m4,compress nogap  s(N r2  F ) ///
          star(* 0.1 ** 0.05 *** 0.01)   b(%6.4f)		  
		  
*面板VAR模型估计
pvar2	  env38 wbgi_rqe deco7 eco13 ,lag(5) soc	 
pvar   wbgi_rqe deco7 eco1  env38,instlags(1/5) fod
est store pvar
esttab pvar ,compress nogap s(N r2  F ) ///
          star(* 0.1 ** 0.05 *** 0.01)   b(%6.3f)
pvarstable,graph
irf create  
pvarirf, impulse(eco1) response(env38)  step(8) mc(50) l(68) byoption(yrescale)
pvarirf, impulse(deco7) response(env38)  step(8) mc(50) l(68) byoption(yrescale)
pvarirf, impulse(wbgi_rqe) response(env38)  step(8) mc(50) l(68) byoption(yrescale)
pvarirf, impulse(env38) response(env38)  step(8) mc(50) l(68) byoption(yrescale)
graph combine rqe.gph eco1.gph eco7.gph env38.gph,row(2)


set scheme s2mono
pvarirf, impulse(wbgi_rqe eco1 deco7 env38) response(env38)  step(10) mc(50) l(68) byoption(yrescale)

